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Abstract 

Deep inelastic scattering data on F% structure function from the various fixed-target ex- 
periments were analyzed in the non-singlet approximation with a next-to-next-to-leading- 
order accuracy. The study of high statistics deep inelastic scattering data provided by 
BCDMS, SLAC, NMC and BFP collaborations was carried out separately for the first 
one and the rest, followed by a combined analysis done as well. For the coupling con- 
stant the following value a s (Af§) = 0.1167 ± 0.0021(totai cxp. error) + | + q QQgg (theor) was 

found, which in this approximation turns out to be slightly less than that obtained at the 
next-to- leading-order, as was generally anticipated. Ditto the theoretical uncertainties 
reduced with respect to those obtained in the case of the next-to-leading-order analysis 
thus confirming earlier observations. 

PACS : 12.38 Aw, Bx, Qk 

Keywords: Deep inelastic scattering; Nucleon structure functions; QCD coupling con- 
stant; NNLO level; 1/Q 2 power corrections. 



1 Introduction 

It goes without saying how it is crucial to know as accurate as possible the parton distri- 
bution functions (PDFs) and the value of the strong coupling constant in order to be able 
to make (relatively) solid predictions for various processes studied in a number of experi- 
ments. Within this realm, the deep inelastic scattering (DIS) of leptons off hadrons serves 
to be a cornerstone process to study PDFs which are universal and feed them further to 
other processes. 

Nowadays the accuracy of data for DIS structure functions (SFs) makes it possible 
to study Q 2 -dependence of logarithmic QCD-inspired corrections and those of power-like 
(non-perturbative) nature in a separate way (see for instance pQ and references therein) 
which is important for the analysis to be performed according to a well defined scheme. 

Until recently a commonly adopted benchmark tool for the analysis happened to be 
there at the next-to-leading-order (NLO) level. However there have already appeared 
papers in which QCD analysis of DIS SFs has been carried out up to the next-to-next-to- 
leading order (NNLO) (see e.g. [2]- [11] and references therein). 

The present paper closely follows the one devoted to the similar study performed at 
NLO level [12] with the major difference in that here we deal with the nonsinglet case 
only because of relative complicacy of the task considered. The singlet part of the anal- 
ysis (combined with the nonsinglet one) will be accomplished in the near future. We 



analyze DIS SF F 2 (x,Q 2 ) with SLAC, NMC, BCDMS and BFP experimental data in- 
volved |13|-|19j at NNLO of massless perturbative QCD. This has become possible thanks 
to the results on both the o? s (Q 2 ) corrections to the splitting functions (the anomalous 
dimensions of Wilson operators) |20] and the corresponding expressions of the complete 
three- loop coefficient functions for the structure functions F 2 and Fl |21j . 

As in our previous paper the function F 2 (x,Q 2 ) is represented as a sum of the leading 
twist F% D { x iQ 2 ) an d the twist four terms Q|: 

F 2 (x,Q 2 ) = Ff CD (x,Q 2 )^+ J ^j . (1) 

While analysing experimental data various corrections must be taken into account. Here 
the nuclear effects, target mass corrections, heavy quark threshold corrections and higher 
twist terms are considered. For details we refer to [121 [22] . 

As is known there are at least two ways to perform QCD analysis over DIS data: the 
first one (see e.g. |23[l24j) deals with Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) 
integro-differential equations [25] and let the data be examined directly, whereas the 
second one involves the SF moments and permits performing an analysis in analytic form 
as opposed to the former option. In this work we take on the way in-between these two 
latter, i.e. analysis is carried out over the moments of SF F 2 (x,Q 2 ) defined as follows H 

M pQCD /twUria/...^ = f 1 x n-2 F pQCD/twist2/...^ q2) ^ ^ 

JO 

and then reconstruct SF for each Q 2 by using Jacobi polynomial expansion method [26] - 
|28j (for further details see |12[ [22] and section [3|) . 



2 A brief theoretical input 

Here we briefly touch on certain aspects of the theoretical part of our analysis. For a bit 
detailed account see [12]. The twist-two DIS SF can be represented as a sum of two terms: 
F^ wist2 (x,Q 2 ) = Fi fS (x,Q 2 ) + F£(x,Q 2 ) , the nonsinglet (NS) and singlet (S) parts. At 
this point let's introduce PDFs, the gluon distribution function fc{x, Q 2 ) and the singlet 
and nonsinglet quark distribution functions fs(x, Q 2 ) and fjvs^, Q 2 ) □: 

/ 

f s (x,Q 2 ) = ^f 9 (x,g 2 ) = F(x,g 2 ) + 5(x,Q 2 ), 

fNs(x,Q 2 ) = u v (x,Q 2 ) - d v (x,Q 2 ) , 

where / is the number of quark flavors (up, down, strange,. . .), V(x, Q 2 ) = u v (x, Q 2 ) + 
d v (x, Q 2 ) is the distribution of valence quarks and S(x, Q 2 ) is a sum of sea parton distri- 
butions set equal to each other. 

There is a direct relation between SF moments ([2]) and those of PDFs 



f N s(n,Q 2 ) = [ dxx n ~ 2 f NS (x,Q 2 ) 
Jo 



For example, in the nonsinglet case it looks [29\: 

M^ S (Q 2 ) = R NS (f) x C^ st2 (n,a s (Q 2 )) x f NS (n,Q 2 ) , (3) 



1 This form was used in [12 [22], too: Eq. (3.32) in [22] should be replaced by Q. 

2 Hereinafter, k — pQG'D, twist denotes the twist two approximation with and without target-mass correc- 
tions (see, for example, [H]). 

3 Unlike the standard case, here PDFs are multiplied by x. 



with 

«,(Q 2 ) = 2^22 (4) 

and C^fg st2 (n^ a s (Q 2 )) are the Wilson coefficient functions. The constant Rns(I) depends 
on the weak and electromagnetic charges and is fixed to be one sixth for / = 4 |29j . 



2.1 Strong coupling constant 



The strong coupling constant is determined from the corresponding solution of the renor- 
malization group equation to an accuracy of C?(10 -5 ) (which is enough for our purposes, 
also we checked that for higher precision the results get no much better). At NLO level 
the latter is given by 
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At NNLO level the strong coupling constant is derived from the following equation: 

a s (Q 2 ) ll + b 1 a s (M 2 ) + b 2 a 2 (M 2 z) 
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for / = 3, 4, 5; A > 0, 



for / = 6; A < 0, 



& are read off from the QCD /3-function: 
= - f3 a 2 - (3ia 3 s - (5 2 a A s + ... 



The equations © and ([6]) allow us to eliminate QCD parameter Aqcd from the anal- 
ysis. However, sometimes it is appropriate to consider it. The coupling constant a s (Q 2 ) 
is expressed through Aqcd (in MS scheme, where Aqcd 

at NLO level 
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A relation between the constant at the normalization point a s (M^) and QCD param- 
eter Aqcd can be obtained from Eqs. ([7]) and ([5D by substituting Q 2 for M|. 
Note that sometimes (see, for example, [2]) the equations 
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are used in the analyses with NLO and NNLO approximations, respectively. These can 
be deduced from the basic equation 
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by expanding an inverse QCD /3-function in RHS of Eq. ([lip (that is l//3(a s )) in powers of 
a s up to 0(a s ) and 0(a 2 ), respectively. The difference between Eqs. Q, ([TP]) and (|7j), © 
can reach 0(1CP 3 ) at Q 2 ~ 1 GeV 2 energies. To avoid uncertainties caused by this 
approach we use in the analyses a numerical solution (with an accuracy of 10~ 5 ) of Eq. flSJ) 
instead. Let's mention in this regard that the approximations given in Eqs. ©, (|10p 
and 0, (jHJ), based on the expansion of inverse powers of In (q 2 /A 2 -^J are very popular 
on the market. They have the following forms: 
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where L^lo = hi(Q 2 : /A^ i0 ) and L = ln(Q 2 /A 2 ) in the NLO and NNLO approximations, 
in order. 
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Figure 1: Difference between a numeric solution to Eqs. and |2P (as precise as O(10~ 5 )) 
and approximate representations given in [W\! and < T73|) /or tne strong coupling constant at NLO 
and NNLO, respectively. 



Note that the difference at NNLO level, shown in Fig. 1, between an approximate ex- 
pression for a s sometimes used in the literature (see, e.g., [30]) and a solution to Eq. (jHJ), 



hence Eq. ([6]), becomes less than that obtained in NLO [31] and still of order 0(1O -3 ) 
(observed in [32j[33]) in the range Q 2 spanning in this analysis (see also discussion in [12]). 
which is comparable with the experimental uncertainties of a s (M|) value obtained from 
the data (see our analyses in Sees. 4 and 5). Although a utilization of the exact tran- 
scendental equation in the NNLO case for deriving the coupling constant appears to be 
not much of a preference over the approximate expression for the latter within an almost 
entire scan range in Q 2 (as opposed to the NLO case), we still prefer to carry out the 
analysis with the former option for it is still an exact (to the order considered) equation 
to be numerically solved to the accuracy desired. 

Also note that the results shown in Fig. 1 were obtained with A(/ = 4) = 200 MeV, 
which is less than those obtained below. Therefore, the actual difference between the 
exact formulas given in Eqs. ([5]), ([6]) and, respectively, their approximations quoted in 
Eqs. (fT2l) . (fT3j) is even a bit more pronounced. 

A starting point of the evolution is taken at relatively large values Qq. There is a 
number of reasons behind that choice, e.g., fewer heavy quark thresholds have to be 
crossed to reach a normalization point, a perturbative approach must be applicable at the 
value of Qq. Besides, impact of higher order corrections derived from PDF normalization 
conditions is the more negligible the higher normalization point is. 



2.2 Q 2 -dependence of SF moments 

The coefficient functions C^ st2 (n, a s (Q 2 Y) is further expressed through the functions 
B J NS (n) which are known exactly |2H I29j 



CffiT(n, a s (Q 2 )) = 1 + a s (Q 2 )B% L s °(n) + aj(Q*)B%$ L " \n) + 0{ai{Q 2 )) 
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The Q 2 -evolution of the PDF moments can be calculated within the framework of 
perturbative QCD (see e.g. [251155]): 
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The function H NS (n, Q 2 ,Ql) up to NNLO may be represented as 
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Here 7^ (n) are the factors before o s in the expansion with respect to the latter of the 
anomalous dimensions ~iNs{n,a s ) (taken in the exact form from |20j). 



4 For the odd n values, the F2 coefficients B 3 NS (n) and Z J NS (n) can be obtained using the analytic continua- 
tion (211 EH ES] . 



2.3 Factorization fip and renormalization /ir scales 



Also, we are set to consider the dependence of results on the factorization fi F and renor- 
malization /iR scales, caused by (see, e.g., pH [38l E5] ) the truncation of a perturbative 
series while doing the calculus. A modification is achieved by replacing a s (defined in 
Eq. (j3J|) in Eqs. (I3|15l) with the expressions in which the scales were accounted in the 
following way: fxp = k F Q 2 , H 2 r = k R ^ F = k R k F Q 2 . 
Then, Eq. (|3|) takes the form: 



M^ S {Q 2 ) = R NS (f) x CTs VI (n,a s (k F Q 2 )) x f NS (n,k F Q 2 ), 



(18) 



and Eq. (|T5j) gets replaced by 

f NS {n,k F Q 2 ) \a s {k F k R Q 
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x H NS (n, k F k R Q' z , k F k R Q$) . (19) 



The functions Cns,H ns are to be obtained from Cns,H ns by modifying the RHS 
of Eqs. (HH H5J) as follows: 
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and in Eq. ([16 
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Actually,while calculating the coefficient functions B^s the renormalization scale depen- 
dence was also taken into account by inserting the term /3q In k R ■ RHS of Eq. (|2ip into the 
expression given in Eq. (|22p and appropriately modifying Eq. (|20p . In this latter case the 
above expressions given in Eqs. (|20p . (j22[) are replaced by the following ones: 
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2.4 Heavy quark thresholds 

Let's now turn to the problem of threshold crossing. We stick to the so-called variable- 
flavor-number scheme, in which any heavy quark of the flavor / is considered to be massless 
and included in the QCD evolution at Q 2 , i.e. Q 2 = Q 2 is the threshold point. The point 



to cross is taken, following (39J 00], to happen at Q 2 ? = m 2 r. H A study into other choices 



5 To be precise, we should have used mj(Q 2 ). However, rrij rather weakly depends on Q 2 in the vicinity of 
Q 2 = m 2 . Thus, hereinafter we adopt m 2 (m 2 ) = mi. The values of the heavy quark masses were taken to be 
those given by the Particle Data Group 2008 42 . 



for threshold crossing and also nowadays popular schemes such as a fixed-flavor one, a 
general-mass variable-flavor-number one and others (see recent paper [H] and discussions 
therein) is deferred to a next paper, with a complete (singlet and nonsinglet) analysis 
carried out. 

Formally, Q 2 evolution does not depend on the specific values of Q 2 ,: a change in 
the initial condition from Q\ 1 to Qq 2 leads only to a change in the normalization from 
Mn (Qo 1) to M^ s (Qq ^)- This property is well reproduced in our analyses (see discus- 
sions at the beginning of Sect. 4). However, the expression for Q 2 evolution does depend 
on the specific values of Qq. 

1. Let Qo be placed in-between the thresholds of / and / + 1 flavors, i.e. Qq = Ql(f)- 
Then, the standard evolution 

M^ s {f,Q 2 ) _ C^ t2 (n,f,af(Q 2 )) f NS (n,f,Q 2 ) 



M n S (f, QlU)) C^ t2 (n, /, al{QD) f NS (n, f, Q 2 (f)) ' 

is correct for Q 2 values between the thresholds of / and / + 1 flavors. Hereafter a( +1 (Q 2 ) 
and a((Q 2 ) denote the coupling constants above and below the threshold Q 2 = Qj + i- 

As it is well-known in the nonsinglet case the coefficient functions beginning at NNLO 
level, and anomalous dimensions starting already with NLO, do depend on the number of 
active quarks. Moreover, starting with NNLO the coupling constant itself is not smooth 
at Qj = m 2 (see 01]). Therefore, we have to deal with the modified equations for 
the latter, which for some heavy quark threshold crossing at Qj+i are found to be of two 
options: 

2. Consider Q 2 evolution above the threshold Q 2 = Qj +1 - Starting from Q 2 = Qj + i, 
it has the above form (|2"3|) with the replacements /—>■/ + 1 and Ql(f) — > Qj + \, that is 

M^ s (f + 1,Q 2 ) _ C^ t2 (n,f + l,af +1 (Q 2 )) ^ i NS (n, f + 1,Q 2 ) _ 



M» S (f + 1, Qj+i) C^ 2 (n, f + 1, d + \Q) +l )) fNs(n, f + 1, Qj +V 

The quantity M^ s (f,Q 2 ) is observable it should be continuous at the threshold: 

M™(f + = M» s (f,Q 2 f+1 ) . (25) 

Therefore, the evolution above the threshold Q 2 = Q 2 + i, in the case of the starting point 
Qq located below the latter, is found to be of the following form: 

M n S U + !>Q 2 ) « st2 (n,/ + W +1 (Q 2 )) ^ f N s(n,f + l,Q 2 ) 
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C^ 2 (n,/,q/(Q 2 +1 )) ^ f NS (n,f,Q 2 +1 ) 
Cftt 2 (nJ,a f s (Q 2 )) X f N s(nJ,Q 2 (f))' 



(26) 



3. Below the threshold Q 2 = Q 2 , we should start from Q 2 = Q 2 and use the expres- 

j+l -> Qj 



sions given in Eqs. (i24"|) and (f26j) with the replacements / + 1 — > f — 1 and Qj + \ — > Q 2 



carried out, i.e. 
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4. By analogy, in the case of two thresholds situated at Q 2 = Q 2 + 2 an< ^ Q 2 = Q'f+i 
and the initial point of the evolution Qq being below the threshold Q 2 = the 
expression is prescribed to be 

M^(/ + 2,Q 2 ) C^ t2 (nJ + 2,4+ 2 (Q 2 )) f NS (n,f + 2,Q 2 ) 
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C^ st2 (n,f,a f s (Q 2 )) X fNs(n,f,Q 2 (f))' 

In the case of two thresholds situated at Q 2 = Q 2 ^\ and Q 2 = Q 2 and the initial point 
of the evolution Qq being above the threshold Q 2 = Q 2 t_i, the rule to follow looks 
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An extension to the case with any number of thresholds is trivial. 
The threshold crossing effect on the coupling constants is implemented according to 
the following equations 
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2.5 Other aspects of the fits 

Analysis's conditions concerning PDF normalization, target mass (TMC) and higher twist 
corrections (HTCs) , as well as nuclear effects remain essentially the same as in our previous 
work [12] so we refer to it for further details, though quoting some salient points. 

The moments fj(n, Q 2 ) at some Qq is a theoretical input to the analysis which is fixed 
as follows. In the fits of data with the cut x > 0.25 imposed only the nonsinglet parton 
density is worked with and the following patametrization at the normalization point is 
used (see, for example, 01!]): 

fjVs(w,Qo) = / dxx n ~ 2 f NS (x,Ql), 
Jo 

fNs(x,Q 2 ) = A NS {Ql){l - x) 6iVs(Q o)(i + d NS {Ql)x) , (31) 
where Ans(Qo)i ^Ars(Qo) an< ^ ^ns(Qo) are some coefficients H. 

6 Here we do not consider the term ~ x a Ns(Q ) j n ^hc normalization of fj^s(x,Qo), because of the cut 
x > 0.25. The correct small- x asymptotics of the nonsinglet distributions is given by Eq. (29) in |12) from the 
corresponding parameters of the valence quark distributions (see Eq. (26) in [T2]) analyzed with allowance for 
the complete singlet and nonsinglet evolutions. 



The distributions of light u and d quarks, f u (x, Qq) = u(x,Qq) and fd(x, Qo) = 
d(x, Qo)i are composed of two components: the valence part — u„(x, Qq) and d v (x, Qq), 
and the sea one — u s (x,Qq) and d s (x, Qq). For the remaining quark and antiquark den- 
sities only the sea parts are retained. Moreover, following [29} 05] an equality of all sea 
parts are assumed with their sum denoted by S(x,Qq). 



3 A fitting procedure 

To cut short this follows along the lines described in the previous paper [12] . Let's here 
just recall salient points of the so-called polynomial expansion method. The latter was 
first proposed in [36] and further developed in [46]. In these papers the method was based 
on the Bernstein polynomials and subsequently used to analyze data at NLO [171 [M] and 
NNLO level [6] [5]. The Jacobi polynomials for that purpose were first proposed and then 
subsequently developed in [261 EH EH] and used in 0-0, [TTlliB]. 

With the QCD expressions for the Mellin moments M^(Q 2 ) analytically calculated 
according to the formula? given above the SF F^x^Q 2 ) is reconstructed by using the 
Jacobi polynomial expansion method: 

Nmax 71 

F 2 k (x,Q 2 )=x a (l-x) b £ Of (x)^4 n) (a,/3)M^. 2 (g 2 ), 

n=0 j=0 

where 0^' 6 are the Jacobi polynomials, a, b are the parameters fitted, and the superscript k 
is defined in the text just before Eq. ([2]). A condition put on the former is the requirement 
of the error minimization while reconstructing the structure functions. 

Since a twist expansion starts to be applicable only above Q 2 ~ 1 GeV 2 the cut Q 2 > 1 
GeV 2 on data is applied throughout. 

MINUIT program [49J is used to minimize two variables 
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where D = dlni^/dlnlnQ 2 . Quality of the fits is characterized by y^/DOF for the 
structure function F%. However, the analysis show that the experimental data for F2 are 
strongly correlated in x and Q 2 ; therefore, it is desirable to have at one's disposal some 
additional characteristics which helps assess the fit quality. From QCD (see subsection 
2.2) it follows that the behaviour F2 ~ (a s (Q 2 )) (see Eq. (|15p . for example) with some d 
values can be taken to be some crude approximation for Q 2 -dependence of the structure 
function. This form is in a sense similar to F 2 ~ (ln(Q 2 /A 2 ))~ d = exp [-din ln(Q 2 /A 2 )] , 
which can be considered as a more appropriate one to use. In other words, the slope 
D ~ —d is approximately Q 2 -independent and, therefore, suffers rather mild correlations 
between x and Q 2 . 

The quantities D th and D exp , corresponding to " experimental data" , can be consructed 
in the following way. Taking several points of the experimental data for F2 with the same 
values of X{, we can parametrize them in the form F2 ~ (ln(Q 2 /A 2 )) di - Xl \ Then, we can 
consider the derivative of this expression with respect to lnln(Q 2 /A 2 ), fitting over each 
subset of Q 2 with the average value of the latter obtained by summing up with the weigth 
l/F(xi, Q 2 ), and then summing over x,. As a result, basically the following expression is 
used 

D = y HQ 2 /A 2 ) dF( Xi ,Q 2 ) 

i F( Xi ,Q 2 ) dln(Q 2 /A 2 )' 

The importance of inclusion of the new characteristics into analysis is shown in Table 4. 
There it is seen that by adding to the fit step-by-step TMC, HTC and systematic errors, its 



quality gradually increases. Indeed, the standard Xg F /DOF demonstrates the fit quality 
improvement, i.e., decreasing from 4.85 to 0.73, while the additional / DOF does it 
even stronger dropping from 55.33 all the way down to 0.71. 

4 Results 

Since there are no gluons in the nonsinglet case the analysis is essentially easier to conduct. 
Hence the cut on Bjorken variable (x > 0.25) imposed where gluon density is believed to 
be negligible. 

We use free normalizations of the data for different experiments. For a reference set, 
the most stable deuterium BCDMS data at the value of the beam initial energy Eq = 200 
GeV is used. With the other data sets taken to be a reference one the variation in the 
results is still negligible. In the case of the fixed normalization for each and all data sets 
the fits tend to yield a little bit worse x 2 ) just as before. 

The starting point of the evolution is taken to be Qq = 90 GeV 2 for BCDMS data as 
well as for overall data and Ql = 20 GeV 2 — for the combined SLAC, NMC and BFP 
data. These Qq values are close to the average values of Q 2 spanning the corresponding 
data. To check for QQ-independence we use also other Qq values: Qq = 2 GeV 2 and Qq 
= 10 GeV 2 . We find that a variation of the results, presented below, is of the order of 
0(1O~ 5 ) for the values of a s (M^) and, therefore, to the accuracy we work in can be said 
to be negligible. 

On grounds of previous knowledge the maximal value of the number of moments to be 
accounted for is N max = 8 [271 12H] (though we check N max dependence like in the NLO 
analysis) and the cut 0.25 < x < 0.8 is imposed everywhere. 

4.1 BCDMS data with carbon, hydrogen and deuterium 
targets 

Analysis commences on with the most precise experimental data [16} \TT\ [T8] obtained by 
BCDMS muon scattering experiment for large Q 2 values. A complete set of data includes 
607 points for the lower cut x > 0.25. As was pointed out earlier the starting point of 
QCD evolution is Qq = 90 GeV 2 . The heavy quark thresholds are taken to be at Q/ = m 2 
(Table 2). An original analysis carried out by BCDMS collaboration (see also |23| ) gave 
(back then) comparatively small values for the strong coupling constant; for example, 
a s (M§) = 0.113 at NLO was quoted in the latter reference. 

Just like in our previous work [T2j an issue with the data systematic errors still remains. 
Let's impose cuts on the kinematic variable Y = (Eq — E)/Eq, where Eq and E are 
lepton's initial and final energies, respectively [50]. Upon excluding a set of data with 
large systematic errors considerably higher values of a s (M§) are obtained and rather 
mild dependence of its values on the choice of Y cut is observed. For more details we refer 

to [12]. 

Impact of experimental systematic errors on the results of QCD analysis as a function 
of Y cut 3, Y cut 4 and Y cut § imposed on data is studied. The following y cuts depending on 
the limits put on x are applied: 



y 


> 0.14 


for 


0.3 < x < 0.4 


y 


> 0.16 


for 


0.4 < x < 0.5 


y 


> Y cu t3 


for 


0.5 < x < 0.6 


y 


> Ycuti 


for 


0.6 < x < 0.7 


y 


> Y cut 5 


for 


0.7 < x < 0.8 



Several cases for the three last conditions, with the cut 0.5 < x < 0.8 imposed on the 
Bjorken variable, are presented in Table 1. 



Table 1. A set ofY cut ^, Y cut ^ and Y cut ^ values used in the analysis 



N Ycut 





1 


2 


3 


4 


5 


6 


Y C ut3 





0.14 


0.16 


0.16 


0.18 


0.22 


0.23 







0.16 


0.18 


0.20 


0.20 


0.23 


0.24 


Y C ut5 





0.20 


0.20 


0.22 


0.22 


0.24 


0.25 



The systematic errors for BCDMS data are given [161 fT71 [T8] as multiplicative fac- 
tors to be applied to F2(x,Q 2 ): fr, fb, fs, fd and fh are the uncertainties caused by the 
spectrometer resolution, beam momentum, calibration, spectrometer magnetic field cali- 
bration, detector inefficiencies and the energy normalization, respectively. 

Each experimental point of the original data set was multiplied by a factor character- 
izing the type of uncertainties under consideration and then the data set modified that 
way was once again fitted along the lines of the procedure given in the previous section. 
The factors f r , fb, f s , fd, fh were read off from [TBI [111 [18]. Absolute differences between 
the a s values for both original and modified data sets are shown in Table 2 in the column 
for a total systematic error estimated in quadrature. There as well given are the number 
of experimental points and a s value for the initial data set. 



Table 2. a s {M%) values for various sets ofY cuts imposed on the data 





number 


X 2 (F 2 )/DOF 


a s (90 GeV 2 ) 


total 


a s (Mf) 




of points 




± stat. error 


syst. error 


± stat. error 





607 


1.06 


0.1523 ± 0.0025 


0.0136 


0.1056 ± 0.0012 


1 


511 


0.96 


0.1671 ± 0.0033 


0.0103 


0.1123 ± 0.0014 


2 


502 


0.96 


0.1680 ± 0.0034 


0.0097 


0.1127 ± 0.0015 


3 


495 


0.95 


0.1685 ± 0.0034 


0.0094 


0.1129 ± 0.0015 


4 


489 


0.95 


0.1701 ± 0.0035 


0.0091 


0.1136 ± 0.0015 


5 


458 


0.94 


0.1719 ± 0.0037 


0.0078 


0.1144 ± 0.0016 


6 


452 


0.93 


0.1729 ± 0.0037 


0.0075 


0.1148 ± 0.0016 



For illustrative purposes let's depict these numbers (to be precise, for a s (M§)) in 
Fig. 2 with NLO results (evaluated in this work) included for comparison. It is seen that 
the value of the coupling is less than in NLO throughout as was generally expected. Also 
note bigger systematic errors with respect to the previous analysis which can presumably 
be ascribed to the scheme of threshold crossing used. 

From the figure one can observe that similar to the analysis done at NLO level the 
values of a s are stable and statistically consistent throughout an entire set of Y cuts 
imposed on the data, though there is a slight trend in the central values of the coupling to 
increase towards higher Ny cut - As in the earlier analysis the case Ny cut = 6 is once again 
most attractive for reducing a total systematic error in ct s by nearly half as much. At the 
same time, increase of the statistical error by 50% is observed just like in NLO case. 

Upon the cuts imposed (in what follows we use the set Ny cut = 6), only 452 points left 
available. Fitting them according to the procedure outlined above the following results 
are obtained: 

a s (90 GeV 2 ) = 0.1729 ± 0.0037 (stat) ± 0.0075 (syst) ± 0.0016 (norm) 

(32) 

a s (M 2 z ) = 0.1148 ± 0.0016 (stat) ± 0.0030 (syst) ± 0.0007 (norm), 
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Figure 2: Effect of systematic errors on the value of the coupling constant for different Y cut 
values in the fits based on the nonsinglet evolution. Data analyzed are BCD MS C 12 , H 2 , D 2 sets 
with the cuts x m i n = 0.25 and those from Table 1 imposed. The starting point for evolution is 
taken to be Qq = 90 GeV 2 . Thresholds of c and b quarks are chosen to be Q 2 = 1.61 GeV 2 
and Qt, = 17.64 GeV 2 , respectively. The inner (outer) error bars show statistical (systematic) 
errors. 

where onwards an abbreviate "norm" denotes the experimental data normalization error 
which comes from the difference of the fits with free and fixed normalizations of BCDMS 
data subsets [T7\ [T8] with different values of the beam energy. Therefore, for the fits of 
BCDMS data in the case of NS evolution and under a condition of minimizing systematic 
errors the following results are obtained: 



where an estimate for the total experimental error comes from the statistical, systematic 
and normalization errors taken in quadrature. Note that this figure is higher than that 
given by BCDMS itself and quoted at the beginning of this subsection. 

Similar to the case of NLO analysis let's scrutinize the dependence of results on a 
maximal number of polynomials N max used in fits. A full set of data is comprised of 452 
points with the Q 2 -evolution starting from Qq=90 GeV 2 . As it can be seen from Table 3 
similar sort of stability of the results still holds in good agreement with [27] . 

As it is seen from Table 3 beginning with N max = 4 the resulting values obtained are 
rather stable; therefore, an average value of the coupling constant can be calculated and 
is found to be a s (M§) = 0.1145. Average deflection is estimated to be 0.0003 and can be 
considered a method error. 



4.2 SLAC, NMC (hydrogen and deuterium), and BFP 
(iron) data sets 



NS evolution analysis is continued with fitting the experimental data obtained by SLAC, 
NMC and BFP collaborations [I3j [TU [151 [19]. A full set of data upon imposing a cut 
x > 0.25 consists of 345 points: 238 SLAC points, 66 NMC points and 41 those of BFP. 
The starting point of the QCD evolution is Qq = 20 GeV 2 and the Q 2 -cut imposed is 



a s (M|) = 0.1148 ±0.0035 (total exp. error) 



(33) 



Q 2 > 1 GeV 2 . For SL AC and NMC data the statistical and systematic errors are combined 
in quadrature. 

Table 3. a s (M§) for various N max values 



Nmax 


X 2 (F 2 )/DOF 


Aslope 

for 6 points 


a, (90 GeV 2 ) 
± 0.0037 


a s (M 2 z ) 
± 0.0012 


3 


1.06 


4.3 


0.1691 


0.1132 


4 


0.96 


5.5 


0.1708 


0.1139 


5 


0.96 


6.6 


0.1702 


0.1137 


6 


0.93 


5.6 


0.1741 


0.1154 


7 


0.93 


4.6 


0.1732 


0.1150 


8 


0.93 


5.0 


0.1729 


0.1148 


9 


0.93 


5.4 


0.1727 


0.1148 


10 


1.05 


7.0 


0.1726 


0.1147 


11 


1.06 


5.4 


0.1716 


0.1143 


12 


1.02 


5.7 


0.1717 


0.1143 


13 


1.10 


5.7 


0.1715 


0.1142 



To illustrate importance of 1/Q 2 corrections the fits of the data are performed in the 
following way. Firstly, one compares the data with the perturbative QCD part of SF F 2 , 
i.e. F2 Wlst2 taken into account. Then, 1/Q 2 corrections beginning with target mass ones 
are added followed by the account for the twist-four terms. As it is easy to read off from 
Table 4 we have unsatisfactory fit when we work with the leading twist part Ft wist2 only. 
Agreement with the data appears to be improving upon including into analysis the target 
mass corrections. Eventually an allowance for the twist-four corrections leads to a very 
good fit of the data. Also, it is seen that the results are considerably spoilt by the neglect 
of systematic errors for SLAC and NMC data, just like those obtained in NLO analysis. 



Table 4. a s (M§) and x 2 f° r various fits with/without TMC, HTC, and systematic 
errors 



N 


TMC 


HTC 


syst. 


X 2 (F 2 )/DOF 


Aslope 


a s (20 GeV 2 ) 


a s {M 2 ) 








error 




for 8 points 


± stat 




1 


No 


No 


Yes 


4.85 


442.6 


0.2260 ± 0.0015 


0.1197 


2 


Yes 


No 


Yes 


2.05 


87.6 


0.2054 ± 0.0014 


0.1139 


3 


Yes 


Yes 


No 


1.47 


14.7 


0.2183 ± 0.0031 


0.1176 


4 


Yes 


Yes 


Yes 


0.73 


5.7 


0.2188 ± 0.0051 


0.1177 



To conclude the following results for x 2 {F 2 ) = 251 and x 2 i ope = ^-7 over 8 points are 
obtained: 

a s (20 GeV 2 ) = 0.2188 ± 0.0051 (stat) ± 0.0084 (syst) ± 0.0025 (norm) 

(34) 

a s (M 2 z ) = 0.1177 ± 0.0014 (stat) ± 0.0035 (syst) ± 0.0008 (norm) . 

The last error ±0.0008 to a s (M z ) comes again from the fits with free and fixed normal- 
izations among different data sets provided by the SLAC, NMC and BFP collaborations. 



Thus, by combining errors in quadrature the fits based on the nonsinglet evolution 
give for the strong coupling constant: 

a s {M 2 z ) = 0.1177 ±0.0039 (total exp. error). (35) 

Looking at the results obtained so far one observes fairly good agreement within errors 
given between the values of a s (M§) derived from the fits of BCDMS data alone and those 
from the fits of combined SLAC, NMC and BFP data. Let's now put all the data together 
and fit them simultaneously. 

4.3 SLAC, BCDMS, NMC and BFP data sets 

Just as above for the BCDMS data the cuts imposed are x > 0.25 along with Y cu t and 
Ny cut = 6 (see Table 1). Then a full set of data consists of 797 points. The starting point 
of the QCD evolution is once again taken to be Ql = 90 GeV 2 . 

Table 5. a s (M^) and x 2 m tie case of the combined analysis 



o 2 


N of 
points 


HTC 


X 2 (F 2 )/DOF 


a s (90 GeV 2 ) ± stat 


a s (Mf) 


1.0 


797 


No 


2.20 


0.1767 ± 0.0008 


0.1164 


2.0 


772 


No 


1.14 


0.1760 ± 0.0007 


0.1162 


3.0 


745 


No 


0.97 


0.1788 ± 0.0008 


0.1173 


4.0 


723 


No 


0.92 


0.1789 ± 0.0009 


0.1174 


5.0 


703 


No 


0.92 


0.1793 ± 0.0010 


0.1176 


6.0 


677 


No 


0.92 


0.1793 ± 0.0012 


0.1176 


7.0 


650 


No 


0.92 


0.1782 ± 0.0015 


0.1171 


8.0 


632 


No 


0.93 


0.1773 ± 0.0018 


0.1167 


9.0 


613 


No 


0.93 


0.1764 ± 0.0022 


0.1163 


10.0 


602 


No 


0.92 


0.1742 ± 0.0023 


0.1154 


1.0 


797 


Yes 


0.98 


0.1772 ± 0.0027 


0.1167 



To verify a range of applicability of perturbative QCD we start with analyzing the 
data without a contribution of twist-four terms (which means Fi = F% ) and perform 
several fits with the cut Q 2 > Q 2 min gradually increased. From Table 5 it is seen that 
unlike the previous analysis [12] quality of the fits starts to appear fairly good already from 
Q 2 = 3 GeV 2 onwards. Except for the order of the approximation at which the analysis 
is performed the basic difference between this analysis and that carried out in [12] is in 
the value of the thresholds that leads to additional increase in the value of the coupling 
constant in comparison to the case of cutting BCDMS data out with large systematic 
errors. Thus, a combination of NNLO approximation and the thresholds taken at Qf = rnt 
rather than Qf = 2my essentially improves agreement between perturbative QCD and 
the experimental data. 

To proceed with comparison, the twist-four corrections are added and the data with 
the usual cut Q 2 > 1 GeV 2 is fitted. It is clearly seen that as in the NLO case here the 
higher twists do sizably improve the quality of the fit, with insignificant discrepancy in 
the values of the coupling constant to be quoted below. 

The following values for the parameters of the parton distribution parametrizations 
for the case corresponding to the last row of Table 5 are obtained: 

A% 2 S = 2.54 ±0.02, A° 2 S = 2.38 ±0.03, A NS = 3.29 ±0.04, A^ e s = 2.35 ±0.17, 
b% 2 s = 4.16 ±0.01, b% 2 s = 4.22 ±0.01, b% s = 4.23 ±0.03, b^ e s = 4.39 ±0.21, 
d% 2 s = 6.08 ±0.17, d° 2 s = 3.89 ±0.12, d NS = 2.02 ±0.19, d N e s = 3.31 ±1.47. 
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Figure 3: Comparison of the HTC parameter h^{x) obtained at LO, NLO and NNLO for 
hydrogen data (the bars stand for statistical errors). 

SLAC NMC BCDMS D 2 data 



0.5 








0.3 



0.4 



0.7 



0.5 0.6 

x 

Figure 4: The same as in Fig. 3 for deuterium data. 
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The parameter values of the twist-four term are presented in Table 6. Note that these 
for H2 and D2 targets are obtained in separate fits by analyzing SLAC, NMC and BCDMS 
data sets taken together. It is seen that the values at NLO and NNLO match within errors 
with an average value being slightly less for the latter. 

For illustrative purposes we visualize them in Figs. 3, 4 where fairly good agreement 
between higher twist corrections obtained at NLO and NNLO is observed, that is in 
agreement with earlier studies (see, for example, |BJ). However, at large x the central 
values of HTCs are a bit decreased at NNLO level. Moreover, in the case of deuterium 
data the HT parameter values in LO for large x are less than those obtained in NLO, in 



contrast to some studies. 



Table 6. Parameter values of the twist- four term in different orders 



X 


LO ± stat 


NLO ± stat 


NNLO ± stat 


hn(x) for H2 


fi4(x) for D2 


hiix) for H2 


hiix) for D2 


hn(x) for H2 


h.4(x) for D2 


0.275 
0.35 
0.45 
0.55 
0.65 
0.75 


-0.266±0.016 
-0.263±0.020 
-0.157±0.037 
0.003 ±0.066 
0.386 ±0.107 
0.901 ±0.187 


-0.273±0.011 
-0.243±0.013 
-0.118±0.022 
0.098 ±0.037 
0.390 ±0.076 
0.549 ±0.122 


-0.258±0.027 
-0.200±0.028 
-0.290±0.023 
-0.316±0.044 
-0.075±0.104 
0.008 ±0.137 


-0.223±0.008 
-0.181±0.007 
-0.005±0.012 
0.244 ±0.020 
0.558 ±0.055 
0.747 ±0.078 


-0.183 ± 0.020 
-0.149 ± 0.028 
-0.182 ± 0.029 
-0.236 ± 0.052 
-0.180 ± 0.135 
-0.177 ± 0.182 


-0.197 ± 0.009 
-0.171 ± 0.015 
-0.033 ± 0.031 
0.142 ± 0.057 
0.295 ± 0.108 
0.303 ± 0.158 



Contrary to [8, 23], we obtain different values of twist-four corrections for the hydrogen 
and deuterium data in NLO and NNLO. Indeed, in the deuterium case, the NLO and 
NNLO corrections have the twist-four corrections insignificantly decreased, whereas in 
the hydrogen case the twist-four corrections are very small at NLO and NNLO. It is quite 
reminiscent of an effect of HTC decreasing in NNLO observed earlier in [3] for F% SF. 



Table 7. Parameter values of the twist-four term in different orders obtained in the 
analysis carried out within a fixed-flavor-number scheme (nf = 4 ) and no cut of BCDMS 
data with large systematics 



X 


LO ± stat 


NLO ± stat 


NNLO ± stat 


h^lx) for H2 


hi(x) for D2 


hi(x) for H2 


h/±(x) for D2 


h^x) for H2 


Ha{x) for D2 


0.275 
0.35 
0.45 
0.55 
0.65 
0.75 


-0.210±0.009 
-0.164±0.010 
0.026±0.016 
0.337 ±0.027 
0.898 ±0.058 
1.508 ±0.113 


-0.193±0.015 
-0.110±0.021 
0.094±0.039 
0.478±0.067 
1.052±0.117 
1.256±0.178 


-0.186±0.010 
-0.149±0.015 
0.009±0.031 
0.260±0.053 
0.719±0.092 
1.179±0.155 


-0.176±0.011 
-0.106±0.015 
0.064±0.030 
0.374 ±0.050 
0.827 ±0.093 
0.918 ±0.134 


-0.163 ± 0.010 
-0.125 ± 0.010 
0.020 ± 0.019 
0.227 ± 0.031 
0.590 ± 0.061 
0.866 ± 0.115 


-0.155 ± 0.012 
-0.087 ± 0.018 
0.066 ± 0.040 
0.324 ± 0.074 
0.667 ± 0.130 
0.606 ± 0.191 



Note that the cut of the BCDMS data, which has increased the a s values (see Fig. 2) 
essentially improves agreement between perturbative QCD and the experimental data. 
Indeed, the HTCs that are nothing else but the difference between the twist-two approxi- 
mation (i.e. pure perturbative QCD contribution) and the experimental data are seen to 
become considerably smaller at NLO and NNLO levels, to compare with NLO HT terms 
obtained in [23 1 and also with the results of analysis obtained within a fixed-flavor-number 
scheme (with a number of flavors fixed to be 4) and no Y-cuts imposed on the BCDMS 
data(see Figs. 5, 6). 

To make it clear with a HTC reduction effect, we perform a few more analyses: 

• within a fixed-flavor-number scheme (FFNS) and rif = 4 (i.e. no thresholds consid- 
ered); 

• no Y-cuts imposed on BCDMS data with large systematic errors (i.e. with Ny cut = 

o); 

• the two above combined. 

As it is seen from Table 7 and Figs. 5 and 6, presented for the last case, without cuts and 
thresholds we reproduce the twist-four corrections obtained in [23J. The corresponding 
values of the coupling constant in NNLO are found to be 

a s {M 2 z ) = 0.1082 for Xsf = 0-96 in the case of # 2 data, 

and 

a s {M 2 z ) = 0.1094 for Xsf = 0-89 in the case of D 2 data. 

It looks like the effect induced by a particular choice of the threshold is small, however 
we plan to study different variants of heavy quark thresholds in our future investigations, 
the type of study carried out in, e.g., [54] . 
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Figure 5: Comparison of the HTC parameter h±(x) obtained at LO, NLO and NNLO for hydro- 
gen data within a fixed-flavor-number scheme (rif = 4) and no Y cuts imposed on the BCDMS 
data (i.e. the case AV ctlt = in Tabl. 1). 
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Figure 6: The same as in Fig. 5 in the case of deuterium data. 



Thus, the combined analysis of SLAC, NMC, BCDMS and BFP data in the case of 
Y-cvXs chosen corresponding to Ny cut = 6 (see Table 1), whenever HTC are not included 
and the Q 2 min cut imposed is 8 GeV 2 (with a free normalization of the data sets), yields 
(for x 2 /DOF = 0.93): 

a s (90 GeV 2 ) = 0.1773 ± 0.0018 (stat) , 



a s (M|) = 0.1167 ±0.0008 (stat) , 



(36) 



and if HTC are included with the cut Q 2 > 1 GeV 2 , correspondingly (x 2 /DOF = 0.98): 

a s (90 GeV 2 ) = 0.1772 ± 0.0027 (stat), 

a s (M z ) = 0.1167 ±0.0010 (stat). (37) 

It is seen that there is no substantial difference between the two, therefore perturbative 
quantum chromodynamics seems to be applicable with the cut Q 2 > 8 GeV 2 imposed as 
this nonsinglet analysis suggests. 

Thus, using the analyses based on the nonsinglet evolution of the SLAC, NMC, 
BCDMS and BFP experimental data for SF Fi with no account for the twist-four correc- 
tions and the cut Q 2 > 8 GeV 2 imposed, we obtain (for x 2 /DOF = 0.93) 

a„(M§) = 0.1167 ± 0.0008 (stat) ± 0.0018 (syst) ± 0.0007 (norm) (38) 

or 

a 5 (M|) = 0.1167 ±0.0021 (total exp.error) . (39) 

Upon including the twist-four corrections, and imposing the cut Q 2 > 1 GeV 2 , the 
following result is found for x 2 /DOF = 0.98: 

a s {M 2 z ) = 0.1167 ±0.0010 (stat) ±0.0020 (syst) ±0.0005 (norm) (40) 

or 

a s (M z ) = 0.1167 ±0.0022 (total exp.error) (41) 

Looking at the results obtained in this section one can note that similar to the NLO 
analysis the central value of the coupling constant a s (M z ) obtained in the fits (NS evolu- 
tion case) of the combined SLAC, BCDMS, NMC and BFP data lie in-between the central 
values of the coupling constant obtained separately in the fits of BCDMS data alone and 
those of SLAC, NMC and BFP data analyzed together. Besides, all the values of a s (M§) 
derived agree within existing statistical errors. 

Within uncertainties, our result for a s (M z ) is also in good agreement with that cited 

in m 

a s (M 2 z ) = 0.1142 ± 0.0023 , (42) 
where a similar analysis of the NS part of the structure function F<i has been performed. 

5 Factorization and renormalization scale depen- 
dence 

In this section the dependence of the results on the different choice of the factorization 
Hf and renormalization /j,r scales are examined. The threshold crossing point is taken to 
be at Q 2 = m 2 because of its substantial role played in the evolution of the coupling con- 
stant [40]. Following the lines of the works [23l [38] we choose just three values (1/2, 1, 2) 
for the coefficients kp and fcjj. 

Results are shown in Table 8. Fits are performed with no account for the higher twist 
corrections, with the number of points equal to 602 (SLAC, BCDMS, NMC, and BFP 
data), with Q 2 min = 8 GeV 2 and a free normalization for different data sets. The change 
in the value of the coupling constant a s {M z ) for various kp and kji values is denoted by 
the difference: 



Aa s (M 2 z ) = a s {M 2 z ) - a s {M 2 z )\ kF=kR=1 



(43) 



Table 8. a s (M z ) for a set of kp and kjt coefficients 



kn 


kp. 


x 2 ^) 


a s (90 GeV 2 ) ± stat 


a s (M|) 


Aa s (M 2 ) 


1 


1 


586 


1773 ± 0018 


0.1167 


o 


1/2 


1 


584 


1734 ± 0017 


0.1150 


-0.0017 


1 


1/2 


585 


0.1717 ± 0.0016 


0.1143 


-0.0024 


1 


2 


600 


0.1845 ± 0.0021 


0.1197 


+0.0030 


2 


1 


592 


0.1829 ± 0.0020 


0.1190 


+0.0023 


1/2 


2 


590 


0.1795 ± 0.0019 


0.1176 


+0.0009 


2 


1/2 


584 


0.1763 ± 0.0018 


0.1163 


-0.0004 


1/2 


1/2 


590 


0.1689 ± 0.0015 


0.1131 


-0.0036 


2 


2 


609 


0.1910 ± 0.0023 


0.1223 


+0.0056 



From Table 8 it follows that the theoretical uncertainties for the maximal and minimal 
values of the coupling constant that correspond to k^ = kp = 2 and Hr = kp = 1/2, 
respectively, are found to be +0.0056 and —0.0036, in order, thus reducing with respect 
to the NLO results obtained earlier [12]. It should be noted that we take into account the 
renormalization scale uncertainty in the expressions for the coefficient functions and the 
respective coupling constants analogously to what was done in [55] . 

Thus, using the analyses with NS evolution of the SLAC, NMC, BCDMS and BFP 
experimental data for SF F2 we obtain for a s (M z ) the following expressions (with no 
account for HTC, Q 2 > 8 GeV 2 and X 2 = 0.93): 

a s {M 2 z ) = 0.1167 ± 0.0008 (stat) ± 0.0018 (syst) ± 0.0007 (norm) 

f +0.0056 ti _. . .... 

+ -0.0036 ( theOT )' ^ 



or 



a s {Ml) = 0.1167 + 0.0021 (total exp.error) + ( (theor). (45) 



6 Conclusions 

In this work the Jacobi polynomial expansion method developed in [26l ETJ [28] was used 
to perform analysis of Q 2 -evolution of DIS structure function F2 by fitting all existing to 
date reliable fixed-target experimental data that satisfy the cut x > 0.25. Based on the 
results of fitting the value of the QCD coupling constant at the normalization point was 
evaluated. Starting with the reanalysis of BCDMS data by cutting off points with large 
systematic errors it was shown that the values of a s (M z ) rise sharply with the cuts on 
systematics imposed. On the other hand the latter do not depend on a certain cut within 
statistical errors. The values a s (M|) obtained in various fits are in agreement with each 
other. An outcome is that quite a similar result for a s (M z ) was obtained in the analysis 
performed over BCDMS data (with the cuts on systematics) and over the data of the rest, 
thus permitting us to fit available data altogether. 

It turns out that for Q 2 > 3 GeV 2 the formulae of pure perturbative QCD (i.e. twist- 
two approximation along with the target mass corrections) are enough to achieve good 
agreement with all the data analyzed. The reference result in NNLO is then found to be 

a s {M 2 z ) = 0.1167 ± 0.0008 (stat) ± 0.0018 (syst) ± 0.0007 (norm), (46) 

Upon adding twist-four corrections, fairly good agreement between QCD (i.e. first two 
coefficients of Wilson expansion) and the data starting already at Q 2 = 1 GeV 2 , where 



the Wilson expansion begins to be applicable, is observed. This way we obtain for the 
coupling constant at Z mass peak at NNLO level: 



a s (M%) = 0.1167 ± 0.0007 (stat) ± 0.0020 (syst) ± 0.0005 (norm) . (47) 

Note that there too is good agreement with the analysis [56] of the combined HI 
and BCDMS data, which was published by HI collaboration. Our result for a s (M z ) is 
also compartible with the world average value for the coupling constant, presented in the 
review [57J 

a s {M 2 z ) = 0.1184 ±0.0007, 
or even more so if it is compared with the recent estimate given by MSTW group |58j : 

a s {M 2 z ) = 0.1171 ± 0.0014 (68%C.L.) ± 0.0034 (90%C.L.) . 

We would also like to note the importance of NNLO corrections in the analyses of DIS 
experimental data. Incorporation of the NNLO corrections have been started already sev- 
eral years ago in various ways. Results are based on the studies of higher order correction 
effects, which can be estimated from the dependence of our results on the factorization 
\ip and renormalization /i/j scales. As was pointed out the values of the theoretical un- 
certainties H, given by this dependence of the results for a s (M z ) are equal to 

a ,,,2m f +0.0056 

A«s(Af z )| theor = | _ 00Q36 . 

For comparison let's quote the analogous numbers obtained at NLO [12J: 

a ,»2m f +0.0070 

Aa s {M z )\ theor = | _ 00Q41 • 

Though the two cases cannot be directly compared, nonetheless some qualitative conclu- 
sions can be drawn. Thus, it is seen that the theoretical uncertainties stay still slightly 
higher than the total experimental error albeit somewhat less than those derived at NLO 
level. Perhaps, this calls for further account of even higher corrections (moreover, maybe 
the ones obtained within approaches different to that we stick with here) and is to be 
given elsewhere. As it was shown in Refs. \55\ 159] . the value of theoretical error should 
decrease approximately by a factor of 2 when the NNLO corrections are accounted for. 
This prediction is hardly observed, which can be attributed to a number of distinctions the 
two analyses bear in part. Though a number of studies, devoted to NNLO QCD analysis 
of the structure functions and appeared in the literature (see [2]- [5], [6j [551 ESI EQ] and 
references therein) in the past, were exploiting back then partially known NNLO QCD 
corrections, it is obvious that in order to analyze experimental data across a whole region 
of x as precise as possible it is necessary to know all NNLO QCD corrections as exact as 
possible. These were evaluated in [2~T] and their exact expressions (rather than the 
approximate expressions given there as well) were used in this paper. 

Concerning the contributions of higher twist corrections in the present work the well- 
known x-shape of the twist-four corrections while going from intermediate to large values 
of the Bjorken variable x is well reproduced. The latter look very similar to those from [23J , 
if no cuts are imposed on BCDMS data with large systematic errors. The latter substan- 
tially reduce the twist-four corrections at NLO and NNLO level, particularly for hydrogen 
data. 



7 It should be mentioned that this analysis was carried out over the data coming from the various experiments 
and in different orders of perturbation theory, i.e., from NLO up to N 3 LO. 

8 As it has already been shown the scale choices \xp — fin — 2Q 2 and hf = = Q 2 /2 give the maximal and 
minimal values of a s (M§) (at the various choices of values kp = 1/2, Uf = 2, ku = 1/2 and ku — 2 separately) 
and thus give main part of theoretical error. 



The next step to take in the study is the consideration of the combined nonsinglet and 
singlet analysis using the DIS experimental data in the full x region and also application 
of some resummation-like Grunberg effective charge method [61] (as it was done in [15] at 
the NLO approximation) and the "frozen" [62] [f| and analytic [63] versions of the strong 
coupling constant (see [631 E3 EZ] for recent studies in this direction) . 

Moreover, we plan to consider also further corrections (i.e. those coming from three 
loops) in the coefficient functions [21j . which permits performing the N 3 LO fits at large 
x values, where the contributions of the corresponding four- loop corrections to the yet 
unknown anomalous dimensions should be negligible. Several N 3 LO fits had already been 
done in [U [HJ [30] . It will be carried out in nearest future with the purpose of studying 
further reduction of theoretical uncertainties. 
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